
#include "datatypes.h"
#include "misclib.h"
#include <limits.h>

#ifndef SERIALQGCM
#include <omp.h> //include openMP
#endif
/* ******************************************************************************************
   define the ``rank'' index struct and force instantiation of required template versions

   The rank class defines an indexing object that is used by multi-dimensional arrays:
   a 1d array is typically indexed by an object of rank<size_t, 1> ---> single integer 
   a 2d array is typically indexed by an object of rank<size_t, 2> ---> vector
   a Nd array would typically be indexed by a rank<size_t, n> --> n-dimensional vector space
 ********************************************************************************************* */

const REAL fNaN = sqrt(-1.0);
const qgcomplex complexNaN( fNaN, fNaN );

// * * force instantiation of commonly needed objects
typedef rank<int,3> triad; 
template class rank<int, 3>;
template class rank<INDEX, 2>;
template class rank<INDEX, 3>;
template class rank<INDEX, 4>;
template class rank<INDEX, 5>;
template class rank<INDEX, 6>;
template class rank<INDEX, 11>;
//-------------------------------------------------------------------------------
template<class T, int s>
rank<T, s>::rank(){
	for(size_t i = 0; i < s; i++) this->list[i] = 0;
};
template triad::rank();
template rank2::rank();
template rank<size_t,2>::rank();
template rank3::rank();
template rank4::rank();
template rank5::rank();
template rank6::rank();
template rank11::rank();
//-------------------------------------------------------------------------------
template<class T, int s>
rank<T, s>::rank(const rank<T,s> & other){
	for(size_t i = 0; i < s; i++) this->list[i] = other.list[i];
};
template triad::rank(const triad&);
template rank2::rank(const rank2&);
template rank3::rank(const rank3&);
template rank4::rank(const rank4&);
template rank5::rank(const rank5&);
template rank6::rank(const rank6&);
template rank11::rank(const rank11&);
//-------------------------------------------------------------------------------
template <>
rank<int, 3>::rank(const int  e1, const int  e2, const int e3){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
};
//-------------------------------------------------------------------------------
template <>
rank<INDEX, 2>::rank(const INDEX  e1, const INDEX  e2){
	this->list[0] = e1;
	this->list[1] = e2;
};
//-------------------------------------------------------------------------------
template <>
rank<size_t, 2>::rank(const size_t  e1, const size_t e2){
	this->list[0] = e1;
	this->list[1] = e2;
};

//-------------------------------------------------------------------------------
template <>
rank<INDEX, 3>::rank(const INDEX  e1, const INDEX  e2, const INDEX e3){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
};
//-------------------------------------------------------------------------------
template<>
rank<INDEX, 4>::rank(const INDEX e1, const INDEX e2, const INDEX e3, const INDEX e4){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
	this->list[3] = e4;
};
//-------------------------------------------------------------------------------
template<>
rank<INDEX, 5>::rank(const INDEX e1, const INDEX e2, const INDEX e3, const INDEX e4, const INDEX e5){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
	this->list[3] = e4;
	this->list[4] = e5;
};
//-------------------------------------------------------------------------------
template<>
rank<INDEX, 6>::rank(const INDEX e1, const INDEX e2, const INDEX e3, const INDEX e4, const INDEX e5, const INDEX e6){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
	this->list[3] = e4;
	this->list[4] = e5;
	this->list[5] = e6;
};
//-------------------------------------------------------------------------------
template<>
rank<INDEX, 11>::rank(const INDEX e1, const INDEX e2, const INDEX e3, const INDEX e4, const INDEX e5, const INDEX e6,
		const INDEX e7, const INDEX e8, const INDEX e9, const INDEX e10, const INDEX e11){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
	this->list[3] = e4;
	this->list[4] = e5;
	this->list[5] = e6;
	this->list[6] = e7;
	this->list[7] = e8;
	this->list[8] = e9;
	this->list[9] = e10;
	this->list[10] = e11;
};
//-------------------------------------------------------------------------------
template<>
rank<INDEX, 11>::rank(const INDEX e1, const INDEX e2, const INDEX e3, const INDEX e4, const INDEX e5, const rank3 & a1, const rank3 & a2){
	this->list[0] = e1;
	this->list[1] = e2;
	this->list[2] = e3;
	this->list[3] = e4;
	this->list[4] = e5;
	this->list[5] = a1[0];
	this->list[6] = a1[1];
	this->list[7] = a1[2];
	this->list[8] = a2[0];
	this->list[9] = a2[1];
	this->list[10] = a2[2];
};
//-------------------------------------------------------------------------------
template<class T, int s>
bool rank<T,s>::operator==(const rank<T,s> & other) const{
	for(size_t i = 0; i < s; i++){
		if( other.list[i] != this->list[i] ) return false;
	}
	return true;
};
template bool triad::operator==(const triad &) const;
template bool rank2::operator==(const rank2 &) const;
template bool rank<size_t,2>::operator==(const rank<size_t,2> &) const;
template bool rank3::operator==(const rank3 &) const;
template bool rank5::operator==(const rank5 &) const;
template bool rank6::operator==(const rank6 &) const;
template bool rank11::operator==(const rank11 &) const;
//-------------------------------------------------------------------------------
template<class T, int s>
T & rank<T,s>::operator[](size_t index){
	dassert(index < s, "Exceeded index max size.");
	return this->list[index];
};
template int & triad::operator[](size_t);
template INDEX & rank2::operator[](size_t);
template size_t & rank<size_t,2>::operator[](size_t);
template INDEX & rank3::operator[](size_t);
template INDEX & rank4::operator[](size_t);
template INDEX & rank5::operator[](size_t);
template INDEX & rank6::operator[](size_t);
template INDEX & rank11::operator[](size_t);
//-------------------------------------------------------------------------------
template<class T, int s>
const T & rank<T,s>::operator[](size_t index) const{
	dassert(index < s, "Exceeded index max size.");
	return this->list[index];
};
template const int & triad::operator[](size_t) const;
template const INDEX & rank2::operator[](size_t) const;
template const size_t & rank<size_t,2>::operator[](size_t) const;
template const INDEX & rank3::operator[](size_t) const;
template const INDEX & rank4::operator[](size_t) const;
template const INDEX & rank5::operator[](size_t) const;
template const INDEX & rank6::operator[](size_t) const;
template const INDEX & rank11::operator[](size_t) const;
//-------------------------------------------------------------------------------
template<class T, int s>
void rank<T,s>::operator=(const rank<T,s> & other){
	if( &other == this ) return;
	for(size_t i = 0; i < s; i++) this->list[i] = other.list[i];
};
template void triad::operator=(const triad &);
template void rank2::operator=(const rank2 &);
template void rank3::operator=(const rank3 &);
template void rank4::operator=(const rank4 &);
template void rank5::operator=(const rank5 &);
template void rank6::operator=(const rank6 &);
template void rank11::operator=(const rank11 &);
//-------------------------------------------------------------------------------
template<class T, int s>
bool rank<T,s>::operator<(const rank<T,s> & other) const{
	for(size_t i = 0; i < s; i++){
		if( this->list[i] < other.list[i] ) return true;
	};
	return false;
};
template bool triad::operator<(const triad &) const;
template bool rank2::operator<(const rank2 &) const;
template bool rank3::operator<(const rank3 &) const;
template bool rank4::operator<(const rank4 &) const;
template bool rank5::operator<(const rank5 &) const;
template bool rank6::operator<(const rank6 &) const;
template bool rank11::operator<(const rank11 &) const;
//-------------------------------------------------------------------------------
template<class T, int s>
bool rank<T,s>::operator>(const rank<T,s> & other) const{
	for(size_t i = 0; i < s; i++){
		if( this->list[i] > other.list[i] ) return true;
	};
	return false;
};
//-------------------------------------------------------------------------------
// Function below is needed for rapid conversion between index sets
// that begin at 1 and those that begin at 0.
template<class T, int s>
void rank<T,s>::operator--(){
	for(size_t i = 0; i < s; i++) (this->list[i])--;

};
template void rank3::operator--();
template void rank11::operator--();
//-------------------------------------------------------------------------------
#ifdef DEBUG // for use with the debugger (gdb)
template<class T, int s>
string rank<T,s>::print() const{
	string rvalue;
	for(size_t i = 0; i < s; i++){
		rvalue += convertint2string( static_cast<int> ( this->list[i] ) );
	}
	return rvalue;
};
template string triad::print() const;
template string rank2::print() const;
template string rank3::print() const;
template string rank4::print() const; 
template string rank5::print() const;
template string rank6::print() const;
template string rank11::print() const;
#endif
/********************************************************************************************
  end struct rank
 ********************************************************************************************/
/* template <class T, int s>
class rankIter {
	private:
		rank<T,s> bounds;
		rank<T,s> value;
	public:
		bool operator==(const rankIter&) const;
		bool operator!=(const rankIter&) const;
		rankIter & operator ++();
		rank<T,s> & operator*();
		bool operator<(const rank<T,s> &) const;
		bool operator>(const rank<T,s> &) const;
		rankIter(const rankIter &);
		rankIter(const rank<T,s> & bounds);
		rankIter(const rank<T,s> & bounds, const rank<T,s> & initialValue);
}; */

// Define rankIter class that allows for single loop iteration of chunks of memory
// that are laid out contiguously according to a set of multidimensional indices 
//-------------------------------------------------------------------------------
template<class T, int s>
rankIter<T,s> & rankIter<T,s>::operator++()
{
	for(int i = s-1; i >= 0; i--){
		if( this->value[i] < this->bounds[i] ){
			this->value[i]++;
			break;
		}
	}
	return *this;
}
template rank4iter & rank4iter::operator++(); 
//-------------------------------------------------------------------------------
template<class T, int s>
rank<T,s> & rankIter<T,s>::operator*()
{
	return this->value;
}
template rank4 & rank4iter::operator*(); 
//-------------------------------------------------------------------------------
template<class T, int s>
bool rankIter<T,s>::operator>(const rank<T,s> & other) const
{
	for(int i = 0; i < s; i++){
		if( this->value[i] > other[i] ) return true;
		if( this->value[i] != other[i] ) return false;
	}
	return false;
}
template bool rank4iter::operator>(const rank4 &) const; 
//-------------------------------------------------------------------------------
template<class T, int s>
bool rankIter<T,s>::operator<(const rank<T,s> & other) const
{
	for(int i = 0; i < s; i++){
		if( other[i] > this->value[i] ) return true;
		if( other[i] != this->value[i] ) return false;
	}
	return false;
}
template bool rank4iter::operator<(const rank4 &) const; 
//-------------------------------------------------------------------------------
template<class T, int s>
rankIter<T,s>::rankIter(const rankIter & other)
{
	this->bounds = (other.bounds);
	this->value = (other.value);
}
// force instantiation
template rank4iter::rankIter(const rank4iter &); 
//-------------------------------------------------------------------------------
template<class T, int s>
rankIter<T,s>::rankIter(const rank<T,s> & bounds)
{
	this->bounds = bounds;
	// rank<T,s> default constructor initializes all elements to 0
}
// force instantiation
template rank4iter::rankIter(const rank4 &); 
//-------------------------------------------------------------------------------
template<class T, int s>
rankIter<T,s>::rankIter(const rank<T,s> & bounds, const rank<T,s> & initialValue)
{
	this->bounds = bounds;
	this->value = initialValue;
}
// force instantiation
template rank4iter::rankIter(const rank4 &, const rank4 &); 
/********************************************************************************************

  define ``multidimarray'' class: 

  multi-dimensional arrays are indexed by an rank<n> index object

Example: multidimarray<float, 2> array2d( 0, rank2(5,5) );
array2d[rank2(1,4)] = 5; // to assign value 5 to matrix element i=2,j=5

class structure:
data is stored in a single 1d array called: T * data;
the array edges define the n-dimensional structure: rank<INDEX, s> edges

 ********************************************************************************************/
template<class T, int s> 
multidimarray<T,s>::multidimarray(){
	this->data = NULL;
	for(size_t i = 0; i < s; i++) this->edges[i] = 0;
};
template multidimarray<qgcomplex,2>::multidimarray();
template multidimarray<qgcomplex,4>::multidimarray();
//-------------------------------------------------------------------------------
template<class T, int s>
const rank<INDEX,s> & multidimarray<T,s>::getedges() const{
	return this->edges;
}
template const rank2 & multidimarray<qgcomplex,2>::getedges() const;
template const rank4 & multidimarray<qgcomplex,4>::getedges() const;
//-------------------------------------------------------------------------------
template<class T, int s>
size_t multidimarray<T,s>::getlinearsize() const{
	size_t size = 1;
	for(int i = 0; i < s; i++){
		size *= static_cast<size_t>( edges[i] );
	}
	return size;
}
template size_t multidimarray<qgcomplex,2>::getlinearsize() const;
template size_t multidimarray<qgcomplex,4>::getlinearsize() const;
//-------------------------------------------------------------------------------
template<class T, int s>
multidimarray<T,s>::multidimarray(const T& initializer, const rank<INDEX, s> & edge){
	this->edges = edge;
	size_t size = this->getlinearsize();
	if(size > 0){
		this->data = new T[size];
		for(size_t i = 0; i < size; i++){
			this->data[i] = initializer;
		}
	}else{
		this->data = NULL;
		for(size_t i = 0; i < s; i++) this->edges[i] = 0;
	}
};

template multidimarray<qgcomplex,2>::multidimarray(const qgcomplex&, const rank<INDEX, 2> &);
template multidimarray<qgcomplex,4>::multidimarray(const qgcomplex&, const rank<INDEX, 4> &);
//-------------------------------------------------------------------------------
template<class T, int s>
multidimarray<T,s>::~multidimarray(){
	delete [] this->data; // data is stored in a single 1d array
}
template multidimarray<qgcomplex,2>::~multidimarray();
template multidimarray<qgcomplex,4>::~multidimarray();
//-------------------------------------------------------------------------------
template<class T, int s>
void multidimarray<T,s>::resize(const rank<INDEX, s> & edge, const T & initializer){
	delete this->data;
	this->edges = edge;
	size_t size = getlinearsize();
	if(size > 0){
		this->data = new T[size];
		for(size_t i = 0; i < size; i++){
			this->data[i] = initializer;
		}
	}else{
		this->data = NULL;
		for(size_t i = 0; i < s; i++) (this->edges)[i] = 0;
	}
};
template void multidimarray<qgcomplex,2>::resize(const rank<INDEX, 2> &, const qgcomplex &);
template void multidimarray<qgcomplex,4>::resize(const rank<INDEX, 4> &, const qgcomplex &);
//-------------------------------------------------------------------------------
template<class T, int s>
const T * multidimarray<T,s>::getptr() const {
	return this->data;
}
template const qgcomplex * multidimarray<qgcomplex,2>::getptr() const;
template const qgcomplex * multidimarray<qgcomplex,4>::getptr() const;
//-------------------------------------------------------------------------------
template<class T, int s>
T * multidimarray<T,s>::getdata()  {
	return this->data;
}
template  qgcomplex * multidimarray<qgcomplex,2>::getdata() ;
template  qgcomplex * multidimarray<qgcomplex,4>::getdata() ;
//-------------------------------------------------------------------------------
template<class T, int s> 
multidimarray<T,s>::multidimarray(const multidimarray<T,s> & other){
	size_t size = other.getlinearsize();
	this->data = new T[size];
	for(size_t i = 0; i < size; i++) this->data[i] = other.data[i];
}
template multidimarray<qgcomplex,2>::multidimarray(const multidimarray<qgcomplex,2> &);
template multidimarray<qgcomplex,4>::multidimarray(const multidimarray<qgcomplex,4> &);
//-------------------------------------------------------------------------------
template<class T, int s> 
multidimarray<T,s> & multidimarray<T,s>::operator=(const multidimarray<T,s> & other){
	size_t size = other.getlinearsize();
	this->data = new T[size];
	for(size_t i = 0; i < size; i++) this->data[i] = other.data[i];
	return *this;
}
template  multidimarray<qgcomplex,2> & multidimarray<qgcomplex,2>::operator=(const multidimarray<qgcomplex,2> &);
template  multidimarray<qgcomplex,4> & multidimarray<qgcomplex,4>::operator=(const multidimarray<qgcomplex,4> &);
//-------------------------------------------------------------------------------
template<class T, int s>
T & multidimarray<T,s>::operator[](size_t i){
#ifdef DEBUG // check for out of bounds

	size_t size = this->getlinearsize();
	for(size_t iter = 0; iter < s; iter++) {
		size*= static_cast<size_t>(edges[iter]);
	}
	dassert(i < size, "Exceeded index max size.");
#endif
	return this->data[i];  
};
template qgcomplex & multidimarray<qgcomplex,2>::operator[](size_t);
template qgcomplex & multidimarray<qgcomplex,4>::operator[](size_t);
//-------------------------------------------------------------------------------
template<class T, int s>
T & multidimarray<T,s>::operator[](const rank<INDEX, s> & pos){
#ifdef DEBUG // assert no index out of bounds
	for(size_t i = 0; i < s; i++){
		dassert( pos[i] < this->edges[i] , "Exceeded edge size.");
	}
#endif
	size_t index1d = pos[0];
	for(size_t i = 1; i < s; i++){
		index1d *= static_cast<size_t>(this->edges[i]);
		index1d += static_cast<size_t>(pos[i]);
	}
	return this->data[index1d];
};
template qgcomplex & multidimarray<qgcomplex,2>::operator[]( const rank<INDEX,2> & );
template qgcomplex & multidimarray<qgcomplex,4>::operator[]( const rank<INDEX,4> & );
//-------------------------------------------------------------------------------
template<class T, int s>
const T & multidimarray<T,s>::operator[](const rank<INDEX, s> & pos) const{
#ifdef DEBUG // assert no index out of bounds
	for(size_t i = 0; i < s; i++){
		dassert( pos[i] < this->edges[i], "Exceeded edge size." );
	}
#endif
	size_t index1d = pos[0];
	for(size_t i = 1; i < s; i++){
		index1d *= static_cast<size_t>(this->edges[i]);
		index1d += static_cast<size_t>(pos[i]);
	}
	return this->data[index1d];
};
template const qgcomplex & multidimarray<qgcomplex,2>::operator[]( const rank<INDEX,2> & ) const;
template const qgcomplex & multidimarray<qgcomplex,4>::operator[]( const rank<INDEX,4> & ) const;
/********************************************************************************************
  end class multidimarray 
 ********************************************************************************************/

/********************************************************************************************
  class regset: 
  container class for model covariance fields, basically a ptr-like class
  that provides for rapid switching between two data sets 1) which stores
  the actual model history, and 2) which is used as temporary storage
  for the model's predictor corrector time integration scheme.
 ********************************************************************************************/
regset::regset(){
	this->data = NULL;
	this->primary_data = NULL;
	this->intermediate_data = NULL;
	this->lsize = 0;
	for(int i = 0; i < 4; i++) this->edges[i] = 0;
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
regset::~regset(){}  // doesn't actually delete anything because the regset is simply a pointer to two multi dimensional arrays (intermediate and primary) 
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
regset & regset::operator=(const regset & od){
	this->data = od.data; 
	this->primary_data = od.primary_data;
	this->intermediate_data = od.intermediate_data;
	this->lsize = od.lsize;
	//#ifdef DEBUG
	this->edges = od.edges;
	//#endif
	return (*this);
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
regset::regset(const regset & od){
	this->data = od.data; 
	this->primary_data = od.primary_data;
	this->intermediate_data = od.intermediate_data;
	this->lsize = od.lsize;
	//#ifdef DEBUG
	this->edges = od.edges;
	//#endif
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::zeromemory(){
	if( ! data ) return; 
	for(int i = 0; i < this->lsize; i++){
		data[i] = 0;
	}
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
bool regset::allocate(const rank4 & edgs){
	this->edges = edgs;
	lsize = 1;
	for(int i = 0; i < 4; i++) lsize *= edges[i];
	this->primary_data = new qgcomplex[lsize]; // new array4d( complexNaN, edges );
	this->intermediate_data = new qgcomplex[lsize]; // new array4d( complexNaN, edges );
	this->setprimary();
	return true;
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex & regset::operator[](const rank4 & pos){
#ifdef DEBUG // assert no index out of bounds
	for(size_t i = 0; i < 4; i++){
		dassert( pos[i] < this->edges[i] , "Exceeded edge size.");
	}
#endif
	size_t index1d = pos[0];
	for(size_t i = 1; i < 4; i++){
		index1d *= static_cast<size_t>(this->edges[i]);
		index1d += static_cast<size_t>(pos[i]);
	}
	return this->data[index1d]; 
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const qgcomplex & regset::operator[](const rank4 & pos) const{
#ifdef DEBUG // assert no index out of bounds
	for(size_t i = 0; i < 4; i++){
		dassert( pos[i] < this->edges[i] , "Exceeded edge size.");
	}
#endif
	size_t index1d = pos[0];
	for(size_t i = 1; i < 4; i++){
		index1d *= static_cast<size_t>(this->edges[i]);
		index1d += static_cast<size_t>(pos[i]);
	}
	return this->data[index1d];
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const qgcomplex & regset::operator[](unsigned int index) const{
	assert( index < lsize );
	return data[index];
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex & regset::operator[](unsigned int index){
	assert( index < lsize );
	return data[index];
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex * regset::get_primary_data(){
	return this->primary_data;
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
int regset::size(){
	return this->lsize; // ->getlinearsize();
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const int regset::size() const {
	return this->lsize; // ->getlinearsize();
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex * regset::get_data(){
	return this->data;
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const qgcomplex * regset::get_data() const {
	return this->data;
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
qgcomplex * regset::get_intermediate_data(){
	return this->intermediate_data;
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::setintermediate() const{
	this->data = this->intermediate_data;
}
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::setprimary() const{
	this->data = this->primary_data;
}
#ifdef DEBUG // functions useful for debugging the program in gdb
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::dump() const{
	this->dumph(getrout());
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::dumpf(const char * fname) const{
	FILE * fout = fopen(fname, "w");
	if(fout){
		this->dumph(fout);
		fclose(fout);
	}else{
		fprintf(stderr, "Could not open '%s' for dumping register.\n", fname);
	}
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FILE * setrout(const char * dev){
	static FILE * rout = NULL;
	if(dev){
		rout = fopen(dev, "w");
	}else if(rout == NULL){
		char * pname = getenv("ROUT");
		if(pname){
			rout = fopen(pname, "w");
		}
		if(!rout) rout = stderr;
	}
	return rout;
};

FILE * getrout(){
	return setrout(NULL);
};

void rclear(){
	fprintf(getrout(), "\033[2J");
};

void rprint(const char * format){
	fprintf(getrout(), "%s", format);
};
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::dumph(FILE * fout) const{
	// dumps the contents of the register to FILE *
	regset temp( *this );
	temp.setprimary();
	fprintf(fout, "DISPLAYING REGISTER\n");
	fprintf(fout, "*this = %p,\tn", this);
	fprintf(fout, "\tDATA = %p,\t", this->data);
	fprintf(fout, "\tPRIMARY_DATA = %p,\t", this->primary_data);
	fprintf(fout, "\tINTERMEDIATE_DATA = %p\n", this->intermediate_data);
	const rank4 & edges = this->edges;
	for(INDEX nki = 0; nki < edges[0]; ++nki){
		for(INDEX ki = 0; ki < edges[1]; ++ki){
			for(INDEX ii = 0; ii < edges[2]; ++ii){
				for(INDEX jj = 0; jj < edges[3]; ++jj){
					const qgcomplex & pdata = temp[rank4(nki,ki,ii,jj)];
					fprintf(fout, "primary [%2i,%2i,%3i,%3i] = (%+12g, %+12g)\n",
							static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(ii), static_cast<int>(jj),
							pdata.real(), pdata.imag()
					       );
				}
			}
		}
	}
	temp.setintermediate();
	for(INDEX nki = 0; nki < edges[0]; ++nki){
		for(INDEX ki = 0; ki < edges[1]; ++ki){
			for(INDEX ii = 0; ii < edges[2]; ++ii){
				for(INDEX jj = 0; jj < edges[3]; ++jj){
					const qgcomplex & pdata = temp[rank4(nki,ki,ii,jj)];
					fprintf(fout, "intermediate [%2i,%2i,%3i,%3i] = (%+12g, %+12g)\n",
							static_cast<int>(nki), static_cast<int>(ki), static_cast<int>(ii), static_cast<int>(jj),
							pdata.real(), pdata.imag()
					       );
				}
			}
		}
	}
	fflush(fout);
};
#endif
// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void regset::deallocate(){
	this->primary_data = NULL;
	this->intermediate_data = NULL;
}

/********************************************************************************************
  end class regset 
 ********************************************************************************************/

/********************************************************************************************
  class twotimevar: container class for model history, also provides interface
  for history lookup
 ********************************************************************************************/
twotimevar::twotimevar(regset & initdata) :
	edges( initdata.edges) 
{
	elements.push_back( initdata ); // add component 0,0
	maxt = 1;
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
//#define get_linear_index(t,tprime) (t*(t+3) >> 1) - (t-tprime)
#define get_linear_index(t,tprime) (t*(t+1) >> 1) + tprime // simplifed the above
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
bool twotimevar::does_element_exist(unsigned int t, unsigned int tprime) const{
	if( tprime <= t ){
		int linear_index = get_linear_index(t,tprime);
		// t*t + 3*t; // formula is (1/2)*(t*t+3*2) - (t-t')
		// linear_index >> 1;
		// linear_index -= (t-tprime);
		return linear_index < this->elements.size();
	}
	return false;
}
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void twotimevar::delvar(unsigned int t, unsigned int tprime){
	int linear_index = get_linear_index(t,tprime); //t*t + 3*t; // formula is (1/2)*(t*t+3*2) - (t-t')
	//   linear_index = linear_index >> 1;
	// linear_index -= (t-tprime);
	rassert( linear_index < this->elements.size() ); // this function is not called frequently
	regset & reg = elements[linear_index];
	delete [] reg.primary_data;
	delete [] reg.intermediate_data;
	reg.deallocate();
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void twotimevar::addeqtime(){
	fatal_stop("Add function body.");
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void twotimevar::addtime(){
	maxt++;;
	for(int i = 0; i < this->maxt; i++){
		regset newset;
		newset.allocate( this->edges ); 
		this->elements.push_back( newset );
	}
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
void twotimevar::addtimes(const vector<int> & tprime_list){
	fatal_stop("Add function body.");
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const regset & twotimevar::getvar(unsigned int t, unsigned int tprime) const{
	int linear_index = get_linear_index(t,tprime); //t*t + 3*t; // formula is (1/2)*(t*t+3*2) - (t-t')
	// linear_index = (linear_index >> 1);
	// linear_index -= (t-tprime);
	return elements[linear_index];
}
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
regset & twotimevar::getvar(unsigned int t, unsigned int tprime){
	int linear_index = get_linear_index(t,tprime); // t*t + 3*t; // formula is (1/2)*(t*t+3*2) - (t-t')
	// linear_index = (linear_index >> 1);
	// linear_index -= (t-tprime);
	return elements[linear_index];
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const regset & twotimevar::getvar_withcgcmexception(unsigned int t, unsigned int t_primed) const{
	// this function should be run in serial
	// It is provided as an interface for the diagnostics scripting language
	//  So, it throws a cgcmexception that can be caught rather than fatally ending the program if a nonexistent set of times
	//  is provided by the user.
	int step = 1; // default advancement
	if( this->does_element_exist( t, t_primed ) ){
		return this->getvar(t, t_primed);
	}
	throw cgcmexception("Could not find requested index set.", -15);
	return elements.front(); 
};
//  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
const int twotimevar::getmaxtime() const {
	return maxt-1;
}
#undef get_linear_index
/********************************************************************************************
  end class twotimevar 
 ********************************************************************************************/

/********************************************************************************************
  class cgcmexception
 *********************************************************************************************/
	cgcmexception::cgcmexception() throw() 
: rvalue(0) 
{}

	cgcmexception::cgcmexception(const string & msg, int r) throw() 
: reason(msg), rvalue(r) 
{}

	cgcmexception::cgcmexception(const char * msg, int r) throw () 
: reason(msg), rvalue(r) 
{};

cgcmexception::~cgcmexception() throw ()
{}

int cgcmexception::id() const throw () {
	return this->rvalue;
};
const char * cgcmexception::what() const throw() {
	return this->reason.c_str();
};

/********************************************************************************************
  end class cgcmexception 
 ********************************************************************************************/


/********************************************************************************************
  class parameter: explicitly defined model parameters are stored in this class 
 ********************************************************************************************/
#define assert_type_check(check_type) if(this->ptype!=check_type) fatal_stop("Invalid parameter access."); 
//-------------------------------------------------------------------------------
bool parameter::IsConstant() const{
	return this->constant;
}
//-------------------------------------------------------------------------------
void parameter::MakeConstant(){
	this->constant = true;
}
//-------------------------------------------------------------------------------
const parameter::type parameter::getType() const{
	return this->ptype;
}
//-------------------------------------------------------------------------------
const vector<REAL> & parameter::getData() const{
	assert_type_check(parameter::Real);
	return ( *(this->values.fvalues) );
};
//-------------------------------------------------------------------------------
void parameter::setValue(const int val)
{
	assert_type_check(parameter::Integer);
	vector<int> & vals = *(this->values.ivalues);
	if( vals.size() > 0 ){
		vals[0] = val;
	}else{
		vals.push_back(val);
	}
}
//-------------------------------------------------------------------------------
void parameter::setValue(const REAL val)
{
	assert_type_check(parameter::Real);
	vector<REAL> & vals = *(this->values.fvalues);
	if( vals.size() > 0 ){
		vals[0] = val;
	}else{
		vals.push_back(val);
	}
}
//-------------------------------------------------------------------------------
void parameter::setValue(const char * val){
	assert_type_check(parameter::String);
	vector<string> & vals = *(this->values.svalues);
	if(vals.size() > 0){
		vals[0].assign(val);
	}else{
		vals.push_back(string(val));
	}
}
//-------------------------------------------------------------------------------
const char * parameter::GetName() const{
	return this->name.c_str();
};
//-------------------------------------------------------------------------------
void parameter::setValue(REAL val, int index){
	assert_type_check(parameter::Real);
	vector<REAL> & vals = *(this->values.fvalues);
	if(index >= 0 && index < vals.size()){
		vals[index] = val;
		return;
	}
	fatal_stop("Attempted to set value of profile out of bounds.");
};
//-------------------------------------------------------------------------------
void parameter::setValue(int val, int index){
	assert_type_check(parameter::Integer);
	vector<int> & vals = *(this->values.ivalues);
	if(index >= 0 && index < vals.size()){
		vals[index] = val;
		return;
	}
	fatal_stop("Attempted to set value of profile out of bounds.");
};
//-------------------------------------------------------------------------------
void parameter::setValue(const char * val, int index){
	assert_type_check(parameter::String);
	vector<string> & vals = *(this->values.svalues);
	if(index >= 0 && index < vals.size()){
		vals[index].assign(val);
		return;
	}
	fatal_stop("Attempted to set value of profile out of bounds.");
};
//-------------------------------------------------------------------------------
void parameter::resize(int nsize) {
	if(this->ptype == parameter::Integer){
		if(nsize <= 0){ 
			this->values.ivalues->clear(); 
		}else{
			this->values.ivalues->resize(nsize, 0);
		}
		return;
	}else if(this->ptype == parameter::Real){
		if(nsize <= 0){ 
			this->values.fvalues->clear();
		}else{
			this->values.fvalues->resize(nsize, 0.0); 
		}
		return;
	}else if(this->ptype == parameter::String){
		if(nsize <= 0){ 
			this->values.svalues->clear();
		}else{
			this->values.svalues->resize(nsize, ""); 
		}
		return;
	}
	fatal_stop("Did not encounter a type.");
};
//-------------------------------------------------------------------------------
const size_t parameter::GetSize() const{
	if(this->ptype == parameter::Integer){
		return this->values.ivalues->size();
	}else if(this->ptype == parameter::Real){
		return this->values.fvalues->size();
	}else if(this->ptype == parameter::String){
		return this->values.svalues->size();
	}
	fatal_stop("Did not encounter a type.");
	return 0;
}
//-------------------------------------------------------------------------------
const int & parameter::getiValue() const{
	assert_type_check(parameter::Integer);
	rassert( this->values.ivalues->size() > 0);
	return ( *(this->values.ivalues) ) [ 0 ];
}
//-------------------------------------------------------------------------------
const REAL & parameter::getfValue() const{
	assert_type_check(parameter::Real);
	rassert( this->values.fvalues->size() > 0);
	return ( *(this->values.fvalues) ) [ 0 ];
}
//-------------------------------------------------------------------------------
const char * parameter::getsValue() const{
	assert_type_check(parameter::String);
	rassert( this->values.svalues->size() > 0);
	return ( *(this->values.svalues) ) [ 0 ].c_str();
}
//-------------------------------------------------------------------------------
const vector<int> & parameter::getiData() const{
	assert_type_check(parameter::Integer);
	return *(this->values.ivalues);
};
//-------------------------------------------------------------------------------
const vector<string> & parameter::getsData() const{
	assert_type_check(parameter::String);
	return *(this->values.svalues);
};
//-------------------------------------------------------------------------------
void parameter::assignData(const vector<int> & idata){
	assert_type_check(parameter::Integer);
	*(this->values.ivalues) = idata;
};
//-------------------------------------------------------------------------------
void parameter::assignData(const vector<REAL> & fdata){
	assert_type_check(parameter::Real);
	*(this->values.fvalues) = fdata;
};
//-------------------------------------------------------------------------------
void parameter::assignData(const vector<string> & sdata){
	assert_type_check(parameter::String);
	*(this->values.svalues) = sdata;
};
//-------------------------------------------------------------------------------
parameter::parameter(){ // default
	this->ptype = Integer;
	this->values.ivalues = new vector<int>;
};
//-------------------------------------------------------------------------------
parameter::parameter(const string & sname, const vector<int> & vals){
	this->constant = false;
	this->name = sname;
	this->ptype = Integer;
	this->values.ivalues = new vector<int>;
	*(this->values.ivalues) = vals;
}
//-------------------------------------------------------------------------------
parameter::parameter(const string & sname, const vector<REAL> & vals){
	this->constant = false;
	this->name = sname;
	this->ptype = Real;
	this->values.fvalues = new vector<REAL>;
	*(this->values.fvalues) = vals;
}
//-------------------------------------------------------------------------------
parameter::parameter(const string & sname, const int & i){
	this->constant = false;
	this->name = sname;
	this->ptype = Integer;
	this->values.ivalues = new vector<int>;
	this->ptype = Integer;
	this->values.ivalues->push_back(i);
}
//-------------------------------------------------------------------------------
parameter::parameter(const string & sname, const REAL & i){
	this->constant = false;
	this->name = sname;
	this->ptype = Real;
	this->values.fvalues = new vector<REAL>;
	this->values.fvalues->push_back(i);
}
//-------------------------------------------------------------------------------
parameter::parameter(const string & sname, const string & s){
	this->constant = false;
	this->name = sname;
	this->ptype = String;
	this->values.svalues = new vector<string>;
	this->values.svalues->push_back(s);
}
//-------------------------------------------------------------------------------
parameter::~parameter(){
	if( this-> ptype == Real ){
		delete this->values.fvalues;
	}else if(this->ptype == Integer ){
		delete this->values.ivalues;
	}else if(this->ptype == String){
		delete this->values.svalues;
	}
};
/********************************************************************************************
  end class parameter
 *********************************************************************************************/
